/*Do-file for "Treating abandoned mine drainage can protect streams cost effectively and benefit vulnerable communities." Run on Stata/SE Version 18.0.
The dataset "Tract_with_AMD_and_EJ.dta" is at the Census tract level, with each row corresponding to a Census tract within Pennsylvania. */

use Tract_with_AMD_and_EJ

*population numbers - note that AMD_quart = 6 are tracts with no streams, which are therefore excluded
total totpop, over(AMD_quart)

*****************
// HH Income
*****************
preserve
collapse (mean) medhhincome (semean) se_medHHI=medhhincome, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuMHHI = medhhincome + 1.96*se_medHHI
gen ylMHHI = medhhincome - 1.96*se_medHHI

twoway (scatter medhhincome AMD_quart if AMD_quart<6) (rcap yuMHHI  ylMHHI AMD_quart if AMD_quart<6)(scatteri 79070 1 "$79,070", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 65608 2 "$65,608", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 62343 3 "$62,343", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 62768 4 "$62,768", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 56116 5 "$56,116", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
 xlabel(none) xtick(1(1)5) xscale(r(0.5 5.5)) xtitle(" ") ytitle("Median Household Income", size(small)) ///
 ylabel(50000(10000)80000) yscale(r(45000 85000))  ylabel(, labsize(small) format(%9.0gc))  legend(off)
graph save Graph "HHIncome_Means.gph", replace

restore
*****************
// Percent Nonwhite
*****************

preserve
collapse (mean) pernonwhite (semean) se_pernonwhite=pernonwhite, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuPNonW = pernonwhite + 1.96*se_pernonwhite
gen ylPNonW = pernonwhite - 1.96*se_pernonwhite

twoway (scatter pernonwhite AMD_quart if AMD_quart<6) (rcap yuPNonW  ylPNonW AMD_quart if AMD_quart<6)(scatteri 20.23 1 "20.2%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 7.48 2 "7.5%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 7.61 3 "7.6%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 9.96 4 "10%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 18.18 5 "18.2%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
 xlabel(none) xtick(1(1)5) xscale(r(0.5 5.5)) xtitle(" ") ///
 ylabel(0(5)20) yscale(r(0 20)) ytitle("Percent Non-White", size(small))  legend(off)
graph save Graph "Race_Means.gph", replace


restore

*****************
// Housing Values
*****************

preserve
collapse (mean) medianhousingvalue (semean) se_avghousingvalue=medianhousingvalue, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuHouseVal = medianhousingvalue + 1.96*se_avghousingvalue
gen ylHouseVal = medianhousingvalue - 1.96*se_avghousingvalue

twoway (scatter medianhousingvalue AMD_quart if AMD_quart<6) (rcap yuHouseVal  ylHouseVal AMD_quart if AMD_quart<6) (scatteri 231041.8258 1 "$231,000", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 157730.8642 2 "$157,700", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 142333.1288 3 "$142,300", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 138792.0245 4 "$138,800", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 117427.6074 5 "$117,400", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
xlabel(1 "No AMD Exposure" 2 "x < 7%" 3 "7%< x {&le} 22%" 4 "22% < x {&le} 51%" 5 "x > 51%", labsize(small) alternate) xscale(r(0.5 5.5)) xtitle(" ") ///
ylabel(110000(40000)230000) yscale(r(100000 240000)) ylabel(, labsize(small) format(%9.0gc)) ytitle("Median Housing Value", size(small)) legend(off)
graph save Graph "HousingValue_Means.gph", replace

restore

*****************
// Percent Renters
*****************

preserve
collapse (mean) perrenters (semean) se_renters=perrenters, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuRenters = perrenters + 1.96*se_renters
gen ylRenters = perrenters - 1.96*se_renters

twoway (scatter perrenters AMD_quart if AMD_quart<6) (rcap yuRenters  ylRenters AMD_quart if AMD_quart<6) (scatteri 28.03287 1 "28%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 20.08141 2 "20%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 21.19397 3 "21.2%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 24.40409 4 "24.4%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 35.16607 5 "35.2", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
xlabel(1 "No AMD Exposure" 2 "x < 7%" 3 "7%< x {&le} 22%" 4 "22% < x {&le} 51%" 5 "x > 51%", labsize(small) alternate) xscale(r(0.5 5.5)) xtitle(" ") ///
ylabel(20(5)40) yscale(r(20 40)) ylabel(, labsize(small) format(%9.0gc)) ytitle("Percent Renters", size(small)) legend(off)
graph save Graph "Renters_Means.gph", replace

restore

// COMBINE THE 4 OUTCOMES INTO 1 GRAPH
graph combine HHIncome_Means.gph Race_Means.gph HousingValue_Means.gph Renters_Means.gph, col(2)
graph save Graph "Means1.gph", replace
graph export "Means1.png", as(png) replace


*****************
// Percent Urban
*****************

preserve
collapse (mean) Pct_urban (semean) se_Pct_urban=Pct_urban, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuPUrban = Pct_urban + 1.96*se_Pct_urban
gen ylPUrban = Pct_urban - 1.96*se_Pct_urban

twoway (scatter Pct_urban AMD_quart if AMD_quart<6) (rcap yuPUrban  ylPUrban AMD_quart if AMD_quart<6)(scatteri 74.40214 1 "74.4%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 40.917071 2 "41%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 46.434556 3 "46.4%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 61.213057 4 "61.2%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 87.526691 5 "87.5%", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
 xlabel(none) xtick(1(1)5) xscale(r(0.5 5.5)) xtitle(" ") ///
 ylabel(40(20)100) yscale(r(35 100)) ytitle("Percent Population in Urban Area", size(small))  legend(off)
graph save Graph "Pct_Urban_Means.gph", replace


restore


*****************
// Weighted Vulnerability Score
*****************

preserve
collapse (mean) Weightedscore (semean) se_Weightedscore=Weightedscore, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuWtdS = Weightedscore + 1.96*se_Weightedscore
gen ylWtdS = Weightedscore - 1.96*se_Weightedscore

twoway (scatter Weightedscore AMD_quart if AMD_quart<6) (rcap yuWtdS  ylWtdS AMD_quart if AMD_quart<6)(scatteri .747 1 "0.8", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 1.3161 2 "1.3", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 1.28208 3 "1.3", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 1.30347 4 "1.3", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 1.23898 5 "1.2", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
 xlabel(none) xscale(r(0.5 5.5)) xtitle(" ") ///
 ylabel(0.75(.25)1.5) yscale(r(0.7 1.5)) ytitle("Energy Transition Vulnerability Score", size(small))  legend(off)
graph save Graph "Vulnerability_Means.gph", replace

restore

*****************
// EJ - EPA with 13 indicators
*****************

preserve
collapse (mean) avg_EJ (semean) se_avg_EJ=avg_EJ, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuEJ = avg_EJ + 1.96*se_avg_EJ
gen ylEJ = avg_EJ - 1.96*se_avg_EJ

twoway (scatter avg_EJ AMD_quart if AMD_quart<6) (rcap yuEJ  ylEJ AMD_quart if AMD_quart<6) (scatteri 5.38 1 "5.4", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 3.38 2 "3.4", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 3.75 3 "3.8", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 4.28 4 "4.3", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 6.03 5 "6", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
xlabel(1 "No AMD Exposure" 2 "x < 7%" 3 "7%< x {&le} 22%" 4 "22% < x {&le} 51%" 5 "x > 51%", labsize(small) alternate) xscale(r(0.5 5.5)) xtitle(" ") ///
ylabel(3(1)6) yscale(r(2.5 6.5)) ytitle("EPA Environmental Justice Score", size(small)) legend(off)
graph save Graph "EJ_Means.gph", replace

restore

*****************
// EJ - CDC
*****************

preserve
collapse (mean) spl_ebm (semean) se_avg_EBM=spl_ebm, by(AMD_quart)
sort AMD_quart
gen x = _n

gen yuEBM = spl_ebm + 1.96*se_avg_EBM
gen ylEBM = spl_ebm - 1.96*se_avg_EBM

twoway (scatter spl_ebm AMD_quart if AMD_quart<6) (rcap yuEBM  ylEBM AMD_quart if AMD_quart<6) (scatteri 7.135 1 "7.1", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 6.396 2 "6.4", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) ///
(scatteri 6.8 3 "6.8", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 7.111 4 "7.1", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)) (scatteri 7.387 5 "7.4", mcolor(edkblue) mlabcolor(black) mlabsize(vsmall)), ///
xlabel(1 "No AMD Exposure" 2 "x < 7%" 3 "7%< x {&le} 22%" 4 "22% < x {&le} 51%" 5 "x > 51%", labsize(small) alternate) xscale(r(0.5 5.5)) xtitle(" ") ///
ylabel(6(0.5)7.5) yscale(r(6 7.5)) ytitle("CDC Environmental Justice Index", size(small)) legend(off)
graph save Graph "EBM_Means.gph", replace

restore


graph combine Pct_Urban_Means.gph Vulnerability_Means.gph EJ_Means.gph EBM_Means.gph,  col(2)
graph save Graph "Means2.gph", replace
graph export "Means2.png", as(png) replace
